Draft version February 2, 2008 

Preprint typeset using 1^1^^ style cmulatcapj v. ^jllj^A, 



TWENTY-ONE CENTIMETER TOMOGRAPHY WITH FOREGROUNDS 

XiAOMiN Wang 

Kavli Institute for Cosmological Physics, Univ. of Chicago, Chicago, IL 60637; xiaomin@cfcp.uchicago.edu and 
Dept. of Physics, Univ. of Pennsylvania, Philadelphia, PA 19104 



Max Tegmark 

Dept. of Physics, MIT, 70 Vassar Street, Rm. 37-626B, Cambridge, MA 02139 and 
Dept. of Physics, Univ. of Pennsylvania, Philadelphia, PA 19104 



o 
o 

(N 
> 

00 

O 

O 
in 
o 

:^ 

o 



X 



Mario G. Santos 

CENTRA, Instituto Superior Tecnico, Universidade Tecnica de Lisboa, Portugal 
Dept. of Physics, Univ. of California at Davis, Davis, CA 95616 

Lloyd Knox 

Dept. of Physics, Univ. of California at Davis, Davis, CA 95616 
Draft version February 2, 2008 

ABSTRACT 

Twenty-one centimeter tomography is emerging as a powerful tool to explore the reionization epoch 
and cosmological parameters, but it will only be as good as our ability to accurately model and 
remove astrophysical foreground contamination. Previous treatments of this problem have focused on 
the angular structure of the signal and foregrounds and what can be achieved with limited spectral 
resolution (channel widths in the 1 MHz range). In this paper we introduce and evaluate a "blind" 
method to extract the multifrequency 21cm signal by taking advantage of the smooth frequency 
structure of the Galactic and extragalactic foregrounds. We find that 21 cm tomography is typically 
limited by foregrounds on scales k <C 1/i/Mpc and hmited by noise on scales k ^ 1/i/Mpc, provided 
that the experimental channel width can be made substantially smaller than 0.1 MHz. Our results 
show that this approach is quite promising even for scenarios with rather extreme contamination 
from point sources and diffuse Galactic emission, which bodes well for upcoming experiments such as 
LOFAR, MWA, PAST, and SKA. 

Subject headings: cosmology: theory — diffuse radiation — methods: analytical 



1. INTRODUCTION 

21 cm tomography is one of the most promising cos- 
mological probes, with the potential to complement and 
perhaps ultimently eclipse the cosmological parameter 
constraints from th e cosmic microwave background 
HRowma.Ti et al\H(M iMi^inn et nimm^. It is also 
a unique probe of the epoch of reionization, which is 
now one of the least understood aspects of modern 
cosmology. There are various techniques to explore the 
epoch of reionization at 5 < z < 20. Apart from the 
CMB (Holder et al. 2003; Knox 2003; Ko gut et a/.ll2003i: 
ISantos et aLi2003.) . radio astronomical measurement of 
21cm radiation from neutral hydrogen has been shown 
theoretically to be a powerful tool to study this period 
(|Madau, Meiksin, & Rocs 1997; Tozzi et al. 2000). Lots 
of work has been done in recent years on various 
theor etical and experimental a spects of 2 1cm radiation 
(e.g. iBarkana k LoeH l2004at iCarilh. Gn cdin. & Owe? 
200a ICarilli et al\ 120041 "^ 'Ciardi & Madau 200E" 
DiMatteo et all I2002t iDimatteo. Ciardi, & Miniap 
20041 'Furlanetto. Sokasian. & Hernauisti | 206T 
Furlanetto^ Zaldarriaga^ fc HernauislT I200I 



I^l< 



I200I; 'Shaver et al' 199^ 'Wvithe & Loe^ l2004albl: 



I'urlanetto^ Zialdarriaga^ & MernauistI lzUU4 

Gnedin fc Shaveil I2003t llliev pJflll I2002all 
Loeb & Zaldarriaga] 12004 : iMcOuinn all 

lOh fc Mackl 120031: 



M orales 12004 
Pen. Wu. fc Petersm] 12004 



ISantos. Coorav. fc Knoxl 



[Zaldarriaga, Furjanetto, & Hernquist 200j 

However, this 21cm tomography technique will only be 
as good as our ability to accurately model and remove 
astrophysical foreground contamination, since the high 
redshift signal one is looking for is quite small and can be 
easily swamped by foreground emission from our galaxy 
or others. With much effort going into upcoming exper- 
iments such as the Mileura Wide-Field Array (MWA)^, 
L0FAR2 (Rottgering 2003), PAST^, and SKA^, aimed at 
gathering redshiftcd 21cm signal from the sky and probe 
the epoch of reionization, it is therefore timely to study 
the foreground problem in detail. 

Although 21cm foregrounds have been discussed 
m some previous papers, (e 0- iDiMatteo et all 120021: 
'Dimatteo. Ciardi. fc Miniati> 120041: iMorales fc Hewit ' 
2003; Oh & Mack 2003; S antos. Coora v. & Knox 200" 
IZaldarriaga. Furlanetto. fc Hernauisti l2003(l . the ques 
tions on how to remove foregrounds and noise from ob- 
servations of 21cm signal, how well it can be done, and 
how reliable it is, are still wide open. Previous papers 
have focused on the angular power spectrum of the sig- 

^ http :/ / web . hay stack . mit . edu /MWA / MWA .html 

^ http://www.lofar.org 

^ http://astrophysics.phys.cmu.edu/ jbp 

* http://www.skatelcscope.org 
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nal, usually assuming a ra ther limited spectral resolu- 
tion tPiMattco et al. 2 001 iDimatteo. Ciardi. fc Miniatl 
|200l, Oh & Mack 2003': 'Santos. Coora v. fc Kno^d 120041: 
I^ldarriaga. Furlanctto. & Hcrnauist 2003)- In this pa- 
per, we develop a method to remove the foregrounds 
along the line of sight, taking advantage of the fact that 
most astrophysical contaminants have much smoother 
frequency spectra than the cosmological signal one is 
looking for. The two approaches are complementary, and 
we will argue that they are best used in combination: 
our technique can be used both to identify point sources 
and other highly contaminated angular regions to be dis- 
carded, and to clean out residual contamination from 
those angular regions that are not discarded. This multi- 
frequency approach is more powerful here than for typical 
CMB applications ( Bennett et al. 2003; Tcgniark' 199Si 
iTegmark. de Qliveir a-Costa. & Hamilton 200^, because 
of the potentially much better spectral resolution, and 
the dramatically oscillating 21cm signal compared with 
smooth foregrounds along frequency direction. 

In this paper, we describe the method for removing 
foregrounds in frequency space, show examples of using 
this method in different scenarios, and discuss its promis- 
ing applications for future experiments. In Section [3 
we introduce the reionization model we use through- 
out the text, then give a brief overview of 21cm emis- 
sion/absorption and computational formalism, on how 
we calculate the 21cm angular power spectrum in Z-space, 
projected ID and 2D power spectra in fc-space, and the 
simulated ID frequency spectrum in real space. In Sec- 
tion 13 we describe our foregrounds removing strategy, 
we also show the foregrounds model we use in our calcu- 
lations. In Section 0] we give several applications of our 
method under different assumptions about foregrounds 
and noise. We summarize our results in Section [S] 

2. REIONIZATION MODEL AND FORMALISM 

The reio nization model we use throughout this pa - 
per is from ^Haiman HoldeHl2f)mtlSa,ntos p.t ffl/.ll200^ 

shown as the solid curve in Fig 3 in ijSantos et al\ 
|2003i^_^Although the most recent results from WMAP 
ijSpergel et QiJ i2006.) favors a lower optical depth, the 
results of this paper are rather insensitive to the de- 
tailed choice of reionization model and associated as- 
sumptions, since we are focused on foregrounds rather 
than the cosmological 21 cm signal. For more informa- 
tion about various reionization models, see, for example, 
(iHaiman fc Holdeil2003lHolder et a/.l2003HSantos et all 

Below we give a brief overview of the 21cm emis- 
sion power spectrum and our calculational method. 
The detailed information on 21cm radiation (emis- 
sion/absorption) phen omena can be found in v arious 
literatures, e.g.. llMadau. Meiksin. fc ReesI 
Santos. Coorav. fc Knoxl l2004t I Shaver et all 



Tozzi et a,l\ 120001 iZaldarriaga. Furlanetto. fc Hernouistl 

mm- 
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3D power spectrum 

The differential antenna temperature observed at 
Earth between the neut ral hydrogen patch a r id the CMB 
can be approximated as lShaver et a/.l l[l999() : lTozzi et al\ 



STf, « 0.016K-(1-|-5)(1 -x) 
h 



0.02 



1-f z 0.3 



10 n„ 



1/2 



(1) 



where S — p/p — 1 is the fluctuation of the density field. 
We write the ionization fraction x as a sum of two terms 
ifSantos et aL.2003^ 

X^Xe{l + S^J, (2) 

where is average ionization fraction and 5^^ is the 
fluctuation of the ionization fraction across the sky. Thus 
the ionization fraction x is not only a function of redshift, 
but also dependent on its position in the sky. 

Assuming (5^1 and Sx^ ^ 1 and neglecting all sec- 
ond and higher order terms, we obtain the 3D power 
spectrum for 21cm emission, 

2 



P3D{k,z) = (0.016K)2 



1 rriuh^ 



V 0.02 

^2g-fe^i?,^.2^ 



H-z 0.3 

10 



X(l - 2Xe +Xi+ h^e xi) Prnatterik) 

where Pmatter{k) is the matter power spectrum, 
the power spectrum for the ionized fraction Pg^ (fc) is 
defined as in ijSantos et a/.ll200^ . 



(3) 
And 



Ps^(k) = b^P„ 



-(k)e 



(4) 



where b is the mean bias weighted by the different halo 
properties . The mean radius R of the ionized patches in 
HII regions is modeled as 

1/3 



R 



1 



1 



R 



(5) 



where Rp is th e comoving size of t he fundamental patch 



Rp ~ lOOKpc 



We assume a cosmological concordance model 
llSeliak et all 12004 iSpergel et all 120031: ITegmark et'dl 
I2004D r»v = 0, flA = 0.71, 17b = 0.047, h = 0.72, 
rig = 0.99, as = 0.9) throughout our calculations. 

2.2. Projected ID power spectra 

The 21 cm signal changes with redshift for two separate 
reasons, one slow and one fast: 

1. The average properties of the Universe {x, Tk, Ts, 
etc.) evolve on a timescale Az ~ 1. 

2. The local properties of the Universe change on 
much smaller scales Az corresponding to the sizes 
and separations between ionized regions. 

Across a very small redshift range zq — Az < z < zo + Az 
where Az < 1, we make the approximation of ignoring 
the former and including only the latter, approximating 
parameters like x, Tk and Tg by their values at zq. This 
enables us to linearizes the relations between frequency 
v, redshift z and comoving radial distance D. 

Making the above-mentioned approximation and ig- 
noring redshift space distortions, the 21 cm signal near 
a given zq has an isotropic 3D power spectrum that 
we can project into ID (radial) power spectra Pm (fc, ^o) 
<Hui. Stebbins. fc Burle 



PlD(fc, Z « Zq) 



1 

2^ 



adial) power sp ectra j 
HT99aiPea,cod<iri99fll) : 

P3c(fc',z«zo)fc'dfc', 



(6) 



3 





k [h/Mpc] 
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Fig. 1. — Line-of-sight ID 21cm power spectra at different red- 
sfiifts and ionization fractions. The thick solid line is power spec- 
trum for neutral IGM at redshift 8.09. 



Fig. 2. — Simulated 21cm signal in frequency space corresponding 
to the z ss 8.09 and x„ « 0.83. 



Figure Q shows the line of sight ID 21cm emission 
power spectra for the fiducial reionization model at differ- 
ent reionization epochs. For comparison, we also plot the 
power spectrum for neutral medium at x — 8.09 (thick 
solid line). 

2.3. Simulated signal in real space from ID power 
spectrum 

We generate and analyze our simulations with fast 
Fourier transforms. The simulated signal in real space 
in the region < x < L is 



1 



JV-l r 

N ^ 



q=0 



A„ 



2'KX 

~ir' 



B„ 



2'KX 

IT' 



where Aq and Bq 
zero mean and standard deviations 



(7) 

are Gaussian random variables with 



AS, = 



v/Pin(fc,zo)/2 = ^JPiD{2'Kq/L,ZQ)/2. N is chosen to be 
a large enough integer that all information from Pioik) 
is included in the summation {kmax — 2'kN/L) . The 
box size L should be small enough so that the range of 
L satisfies Az <C 1. Figure |21 shows our simulated 21cm 
signal versus frequency around 155 MHz, corresponding 
to an epoch around redshift zq — 8.09. 

If we plotted the observed signal 
evant frequency range, the expected 
from foregrounds would lie far above 
logical signal shown in. The key to 
cosmology is therefore removing foregrounds using 
multi-frequency information, as emphasize d by, e.g., 
(jZaldarriaea. Furlanetto. fc Hernauisq I2003D . and we 
now turn to this subject. 



in the rel- 
contribution 
the cosmo- 
doing 21cm 



METHOD FOR FOREGROUND REMOVAL IN 
FREQUENCY SPACE 



Because of its small frequency cross-correlations, the 
21cm signal is oscillating dramatically along the fre- 
quency direction. The foregrounds, on the other hand, 
are generally quite smooth over the short frequency range 
we consider. This slowly varying nature of the fore- 
grounds comp ared to the signal is a grea t advantage when 
removing it l|( j}nedin fc Shaved 1 200.'^ iMcQuinn et a,l\ 



120051 iMorales. Bowman, fc Hewitij I2005D . and it is the 

main reason that our foreground removal method works 
so well. 

Our method described here is insensitive to the reion- 
ization model and the redshift range we choose, since we 
are focused on foregrounds rather than the cosmological 
21 cm signal. 

3.1. Foreground removal method 

Our basic approach is to subtract foregrounds sepa- 
rately in each angular direction in the sky, by fitting 
their total intensity dependence of frequency by a log- 
log polynomial. Note that since we are fitting the total 
foreground spectrum separately pixel by pixel (fitting not 
only for the amplitude but also for the spectral index and 
the running of the spectral index), we are unaffected by 
the possible complication of huge variations of the fore- 
ground spectral index across the sky. (If the foregrounds 
would lack both frequency coherence and spatial coher- 
ence, i.e., fluctuate randomly with both frequency and 
position, then we would be unable to identify and re- 
move them and could merely average them down like we 
do with detector noise.) 

There are two separate steps in our analysis: 

1. Simulation 

2. Cleaning 

We treat them as completely independent. In other 
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words, our cleaning algorithm is blind, containing no in- 
formation about the foreground and noise model used in 
the simulation step. It is entirely specified by the single 
integer m giving the order of the log-log fitting polyno- 
mial. 

In the simulation step, we simulate for each pixel the 
total observed frequency spectrum yi at n different log- 
frequencies Xi = \%{vi), i — l,...,7i. This simulated to- 
tal signal y, = Ig + Ily„ + P^J, in- 
cludes 21cm signal /licmi synchrotron emission fore- 
ground free- free emission foreground J^^, point 

source foreground Ip^ and detector noise I^^^. We test a 
variety of different assumptions for foregrounds and noise 
in this step. 

Then we group the y^'s into an n-dimensional vector 
y, and group the x^'s and their powers into an n x to 
matrix X. So that the data can be modeled as 



y = Xa + n, 



(8) 



where the m-dimensional parameter vector a 
parametrizes the foreground contributions. It is 
what we need to find out in the cleaning step. In 
equation (jSJ, n is the part left in the total signal that 
can not be fitted by the parameters in a, including the 
contribution from signal, detector noise and residual 
foregrounds. 

In all our calculations throughout this paper, we take 
this fitting polynomial to be quadratic, i.e., fit the to- 
tal foregrounds as a single running power law. Equiv- 
alently we fit the log intensity of the foreground as 
Ig / = a^ + a2^gi^ + ai (Ig i^)^ • That is to say, the number 
of fitting parameter in our computation is always m = 3. 
We found that this improved noticeably over m = 2, 
whereas increasing to m = 4 gave essentially no further 
improvement. 

In the cleaning step, we find our best fit parameter 
vector a b y minimizing = (y ~ Xa)*N^^(y — Xa), 
obtaining (jTegmarkll 19971) 



a= [X*N-iX]-iX*N-V, 



(9) 



where N is the covariance matrix of the contributions 
from the detector noise and 21cm signal. We then sub- 
tract the total fitted foreground contribution Xa from 
the simulated measurement vector, thus obtaining what 
we will refer to as the cleaned signal y. 

Although this cleaning technique is only optimal if N 
is known and the c ontributions fro m noise and 21cm sig- 
nals are Gaussian l)Tegmarklll997|) . we use equation ^ 
anyway and quantify the residual noise using our sim- 
ulations. Since equation ^ minimizes the rms residua l 
even in the presence of non-Gaussianity ljTegmarklll997lj) . 
it is a robust general-purpose fit that does not require 
detailed foreground or signal modeling. We simply set 
N = I, the identity matrix, which will be essentially opti- 
mal if white detector noise dominates. If desired, this can 
be further improved by modeling the foreground power 
spectrum found in real data and iterating. 

Since the cleaning step uses a single polynomial in log- 
log, it cannot fit exactly a simulation including detec- 
tor noise or more than one foreground component (since 
adding the exponential of two polynomials does not give 
the exponential of a polynomial). We will see that this 
simple cleaning algorithm is nonetheless very successful. 



able to fit any of our foreground models well over the 
limited frequency range that is relevant. 

3.2. Foreground and detector noise models 

The foregrounds we consider in this paper include 
Galactic synchrotron emission, Galactic free-free (ther- 
mal) emission, and extragalactic point sources. For 
more information on models of fore grounds in the 
100 MHz range, please see , e.q., iPlMatteo et al 



1982; 



iUU MMz range, please see, e.g.. iiJJiiviatteo et at. 
2002t iDimatteo. Ciardi. fc Miniatil 120041: iHaslam al 



iHaverkormKatgert"^ de Bruvn" 
Morales fc Hewitt 2003: Oh fc Mack 



2003; 

, 2003; 

Platania et al\ 120031 ISantos. Coorav. fc Kno:?d 12004 ; 
Shaver et all Il999j) . As emphasized above, the fore- 
ground models we describe here are used only in our 
simulation step, not for the cleaning process. 



3.2.1. Galactic synchrotron radiation 

For Galactic synchrotron emission, which probably 
causes most contamin ation of all foregrounds ( perhaps of 
order 70% at 150MHz ijPlatania e t al. 2003; Sh aver et all 
11999(1 ). we assume its intensity to be a running power law 
in frequency 



^syn 



(10) 



with a spectral index asyn = 2.8 ijTegmark et a,l\ 
12000^ and a spectral index "running" Aa^un = 0.1 
iHay erkorn. Katgcrt, fc de Bruvn 2003: Platani a et all 
l2003t Igha^er et all ll99Et iTegmark et all l2000|) . Here 
i/^ = 150MHz. We assume the amplitude of the syn- 
chrotron foregr ound to be ^.si/n 



335.4K, an extrap- 



olation from fHa verkorn. Katgert. fc de Bruvn| 12003'; 
l^gmark et al. 200(|). We also explore other normaliza- 
tions Asyn orders of magnitude higher than the value we 
define here in our calculations. Similarly we try other 
values of spectral index and spectral running index in 
the calculation to test the robustness of our method. We 
discuss the details in Section^ 

3.2.2. Galactic free-free emission 

We model the Galactic free-free emission (which 
might contribute a con tamination of ord er 1% con- 
tamination at 150MHz fShaver [1^ 991) ) as a run- 
ning power law as we ll (Havcrkorn^ Katgert^ fc de Bruvi] 
l2003tlPlatania7nDl2003HTegmaTk et all\2mf. 



Iff — Aff I — 



-Off-Aoff Ig 



(11) 



where a// = 2.15, Aaff = 0.01, and Aff = 33. 5K, ex- 
trapolated from (Haverkorn. Katgert. fc de Bruvn..200^i 
ITegmark et'!Ii\\2mCl^ . 

3.2.3. Extragalactic point sources 

Point sources have been estima ted to cause about 30% 
of the contamination at 150MHz l|Shaver et a/.ll99^ and 

are typically less smooth in frequency than the Galactic 
foregrounds. When looking in a given direction, we are 
observing the same point sources as we change frequency, 
so there are not small-scale fluctuations in the same sense 
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as when we change observing directions"' . A serious com- 
phcation compared to the Galactic synchrotron and free- 
free cases is that when we observe many point sources in 
a pixel, they can each have quite different spectral in- 
dices, possibly making their combined intensity a quite 
complicated function of frequency. 

One approach would be to model this complicated 
function as a running power law over the narrow fre- 
quency range involved, just as we did for the synchrotron 
and free-free foregrounds: 

/ C \ P / ,, \ -apa-AOpa Ig ■" 



where a„ 



\mjy 
2.81 



(12) 



Aap, = 0.25 and /3 = 0.125 

ijTegmark et a/.ll2000j) . 

However, we wish to be as conservative as possible in 
this paper, and therefore adopt a more complicated point 
source model in our simulations. Instead, we therefore 
simulate a large number of random point sources i — 
1, ... in the pixel that we are considering and sum their 
intensity contributions in Kelvin: 



where 

(Ips) 



dB 

It 



150 MHz\"' 



dB 



,dN 
"dS 



dS 



150 



'{Ips), (13) 



f{a) da 



. (14) 

is the average value of point source foreground inten- 
sity. Conversion factor dB/dT = 6.9 x lO^mJy/K. The 
assumed sky area per pixel is approximately fisky = 
12 arcmin^. In equation ifT^ . S* is the flux of the i"^ 
point source at 150 MHz. It is generated randomly from 
the source count distribution dN/ dS = 4(S'/lJy)^^ '''^ 
l|Dimatteo. Ciardi. fc Miniatill200^ , truncated at a max- 
imum flux S'max = Scut = O.lmJy, above which we as- 
sume that point sources can be detected and their pix- 
els discarded. In other words, when we talk about the 
contamination from point sources below, we refer only 
to the contribution from unresolved point sources. To 
avoid having to generate infinitely many point sources, 
we also truncated the distribution at a minimum flux 
'5'min = lO^'^mJy, since we find that the total flux con- 
tribution has converged by then. We generate a^, the 
spectral index of the z*'' point source, randomly from the 
Gaussian distribution 



/(«) 



1 



■ exp 



{a - ao) 



2al 



with spectral index a in the range of [ag — 
Aa], where Aa — ha a- To be conservative 
the spectral index vary in a fairly large region, 
through our calcuations. 



(15) 

Aa,ao + 
we allow 
10, 



^ There is, however, the subtle effect of off-beam point 
sources dimming towards higher frequencies because tfie beam gets 
narro wer (^h & Mack 2003; Zaldarriaea. Furlanctto. & Hernauist 
12003) . For the narrow frequency intervals A In that we are con- 
sidering, this effect will be around the percent level for individual 
sources, in the same ballpark as the intensity change due to the fre- 
quency dependence of the emission mechanism. This means that 
it will not imprint sharp spectral features in the total foreground 
emission, and should be well fit by our blind cleaning algorithm. 
We have not included this complication in the present analysis — 
it would be worth incorporating it in a more detailed foreground 
analysis, particularly in one including explicit modeling of the sky 
pixelization. 



3.2.4. Detector noise 

We treat detector noise as white noise. In the Rayleigh- 
Jeans limit, the rms detector noise in a pixel can be ap- 
proximated as 



A2 



-B 



A2 S 



(16) 



2kB 2kB A 

where kg is Boltzmann constant, A is redshifted wave- 
length of 21cm emission. The specific brightness B is 
related to point source sensitivity S by dividing it with 
pixel area A. 

At redshift 8.47, v = 150MHz, A = 2m, with LOFAR 
virtual core configuration^, for a 5.2' pixel with 4MHz 
bandpass and one hour integration, sensitivity S is ap- 
proximately 0.17mJy, from equation H16|) . we get 



= 108mK (IM^) 



0.5 



1 hour\ 
t j 



0.5 



(17) 



where Ai^ is channel width and t is total integration time. 
Similarly for MWA experiment^, a 4.6' pixel with 32MHz 
bandpass and one hour integration, point source sensitiv- 
ity S — 0.27mJy, so we get MWA detector noise 



,MWA 



= 218mK 



/32MHz\ 



0.5 



1 hour^ 



0.5 



t 



(18) 



We should mention that although at 4MHz bandwidth, 
the sensitivity for MWA is worse than LOFAR, MWA has 
a larger bandpass and field of view. This larger field of 
view leads to vastly more pixels, which is an advantage 
for foregrounds removal, as we will see in later sections. 
The detector thermal noise is only one of the many con- 
cerns in the experiment, such as calibration, systematics, 
etc. Thus it should not be considered as the only criterion 
to judge an experiment. 

The 1-dimensional power spectrum of the detector 
noise can then be written as 



Pdet = 2Traj^. 



(19) 



In our simulation, we consider two scenarios. One sce- 
nario assumes a fiducial future experiment with Gaus- 
sian random detector noise down to ax = ImK level. 
The other scenario assumes a currently achievable de- 
tector noise level of ~ 200mK. This is based on equa- 
tion |(T7|) and equation UHl) for the LOFAR and MWA ex- 
periments, assuming 1000 hours of integration time and 
4kHz ~ 8kHz frequency resolution, respectively. 

4. RESULTS 

As we showed previously in FigureEl the signal wiggles 
rapidly with frequency. This is the key advantage of re- 
moving foregrounds in frequency space, since foregrounds 
are typically relatively smooth functions of frequency. 

We simulate the 21 cm signal as a Gaussian random 
field, although in reality, the signal is of course highly 
non-Gaussian. We make this Gaussianity approximation 
for simplicity, since the key quantity that we are inter- 
ested in (the power spectrum of the residual noise and 
foregrounds) depends mainly on the power spectrum of 
the signal, foregrounds and noise, not on whether the 
statistics are Gaussian or not. 

^ http://www.lofar.org 

^ http :/ / web . hay stack . mit . edu /MWA / MWA .html 
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4.1. Baseline example 1 — long term potential 
(noised signal) 

The results for the basehne example with noise much 
smaller than the signal are shown in Figure |31 The top 
panel shows the total contaminant in a pixel including 
Galactic synchrotron radiation, Galactic free-free emis- 
sion, extragalactic point sources and detector noise with 
a = ImK, which is the fiducial value for a future gen- 
eration experiment. The foregrounds are modeled as in 
the previous section with parameters (given in the figure 
caption) corresponding to rather pessimistic assumption 
about the foreground properties. 

The total foreground is so huge that the top panel looks 
really similar to the middle panel (which includes the 21 
cm signal) . From the middle panel, we can not really tell 
if there is any signature of 21cm emission. However, the 
bottom panel shows that the foregrounds can be effec- 
tively removed, with the residual (recovered 21 cm sig- 
nal minus the input "true" signal) being more than five 
orders of magnitude below the foregrounds in amplitude. 
By transforming the detector noise and the residual from 




Fig. 3. — Spectrum in a single pixel before and after foreground 
cleaning. The top panel shows total contaminant signal, consisting 
of synchrotron radiation (Asyn = 335. 4K, asyn = 2.8, Aasy„ = 
0.1), free-free emission foreground {Aff = 33. 5K, ajj = 2.15, 
A«j^j^ = 0.01), extragalactic point sources {aa = 10) and detector 
noise (cr = ImK). The middle panel has the cosmological 21cm 
signal added. The bottom panel shows the recovered 21cm sig- 
nal (blue curve) compared with the true simulated signal (black 
curve) and the residual (recovered minus simulated 21cm signal; 
red curve). The three horizontal black dashed lines correspond to 
-0.004K, OK, and 0.004K, respectively. (Note the different vertical 
axis limits.) The small-scale wiggles in the residual are detector 
noise, whereas the smoothed parabola-shaped component of the 
residual is the error in the foreground fitting. 

the bottom panel of Figure|3|back to A:-space, we are able 
to compare them with the 21cm ID power spectrum, 
shown in Figure ^ Before foreground cleaning, the to- 
tal contaminant (blue dotted curve) is seen to dominate 
over the 21cm power spectrum (black solid curve). The 



foreground power spectrum is seen to rise towards the 
left, reflecting its rather smooth frequency dependence. 
After foreground cleaning, the residual contaminant (red 
dashed curve) is significantly below both the original con- 
taminant except on scales k <C 0.1/i/Mpc. The flat sec- 
tion of the residual power spectrum for k > 1/i/Mpc is 
seen to correspond to detector noise. 

The three vertical lines in Figure 0] correspond to dif- 
ferent minimum channel width for different upcoming ex- 
periments. From left to right, they are O.lMHz (fiducial), 
8kHz (for MWA«) and 4kHz (for LOFAR^). Information 
to the right of this minimum channel width line is lost, 
roughly corresponding an exponential blow-up of the ef- 
fective detector noise. In other words, the effective detec- 
tor noise goes much higher above the signal to the right 
of those lines so that little information can be extracted 
there. So in order to take advantage of our method of 
foreground removal, the channel width needs to be small 
enough to reach the noise-dominated region. 




k [h/Mpc] 



Fig. 4. — ID power spectra of 21cm signal (solid black curve), 
total contaminant from Figure|21 (dotted blue curve), residual con- 
taminant after foreground cleaning (dashed red curve) and detec- 
tor noise alone (dotted black curve). The horizontal solid blue 
line is the white noise power spectrum used for the detector noise 
simulation. The three vertical lines correspond to the different 
channel width for different experiments. From left to right, they 
are O.lMHz (fiducial, short dash - long dash blue line), 8kHz (for 
MWA, dot - long dash black line), and 4kHz (for LOFAR dot - 
short dash red line). 

In other words, 21 cm tomography is limited mainly by 
foregrounds for k <^ 1/i/Mpc and limited mainly by noise 
for k 3> 1/i/Mpc. To take full advantage of their sensitiv- 
ity by pushing residual foregrounds down to the detector 
noise levels, experiments should therefore be designed to 
have a channel width substantially smaller than 0.1 MHz. 
Such small channel widths are realistic and achievable for 

* http://web.haystack.mit.edu/MWA/MWA.html 
^ http://www.lofar.org 
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upcoming experiments. Since the analysis can now prac- 
tically be done by dedicated high speed electronics, even 
if the software solution was not fast enough. 

To test the robustness of our foreground cleaning tech- 
nique, we repeated the above analysis for a wide range 
of foreground models with the same noise level. 

First we tested a suite of models with only detector 
noise and synchrotron radiation, changing the values of 
the parameters defined in equation IjlUI) . Most of the re- 
sults were similar to those shown in Figure |21 Increasing 
the synchrotron amplitude parameter Asyn all the way 
up to completely unrealistic value lO^K had essentially 
no effect, because in this case the simulated spectrum had 
the exact same shape as the model fit for in the cleaning 
step. Likewise, changing the spectral index asyn over the 
extreme range —80 < asyn < 20 had little effect. Com- 
plicating the synchrotron spectrum with a running of the 
running term Aa^^„ lg^(-^), so that the intensity of the 
synchrotron foreground can be written as 

/ ^ \ -asy,^-^asy,^ Ig ~ Aq Ig^ ( ) 

■^syn — ^syn \ 1 i 

(20) 

still caused a negligible increase of the residual as long 
as lAofsJ^jJ < 50. This is a reflection of the fact that 
over a fairly narrow frequency range, even a more com- 
plicated function can be accurately approximated by a 
parabola in log-log, i.e., by the first three terms of its 
Taylor expansion. 

Making variations around the baseline case of multiple 
foregrounds, similarly we tried numerous examples with 
either higher foreground amplitudes, different spectral 
indices or larger running of the spectral indices, again 
obtaining results similar to the ones shown in Figure |3I 
For example, for the point source foreground, we tried 
different values for parameter (Tq in equation IjlSI) from 
(Tq — 0.2 up to (Ja = 100. We found that as long as 
(Tq < 60 (a conservative cut), the residuals are rather 
insensitive to the distribution of point source spectral 
indices. 

All our tests show that as with astrophysically plausi- 
ble foreground amplitudes, the effectiveness of our sim- 
ple "blind" cleaning method is almost independent of the 
number, shape and amplitude of relatively smooth fore- 
grounds. 

The basic reason for this robustness is easy to under- 
stand. As long as the total foreground signal can be well- 
approximated by our fitting function (a log-log parabola) 
over the small frequency interval in question, then the 
main contribution to the residual will not be the fore- 
grounds, but the amplitude of this log-log parabola con- 
tributed by random detector noise. Our numerical exam- 
ples show that the residual indeed does have roughly the 
shape of our fitting function, not the shape of the main 
leading order contribution from residual foregrounds (the 
next term in their Taylor expansion, i.e., a cubic term). 

4.2. Baseline example 2 — near term situation 
(noised signal) 

The detector white noise we assumed in previous ex- 
ample is small comparing to the signal, in which case we 
can subtract the foreground easily for each pixel. Nev- 
ertheless, that level of noise might not be achieved un- 
til future next generation experiments. For upcoming 



experiments, as we showed in equation H17|l and equa- 
tion H18|l . the detector noise is far above the 21cm signal. 

In this section, we study the close term case by as- 
suming detector noise a = 200mK, 200 times larger than 
previous example, using the same baseline foreground 
model as in previous section. 

One would expect that for each individual pixel, the 
huge white noise destroys much of the information about 
the 21 cm signal and also obscures the frequency depen- 
dence of the foregrounds, makeing them harder to fit and 
remove, and that the residual will be noise-dominated. In 
this scenario, single pixel cleaning is not enough for clean- 
ing purposes and multiple pixels are needed to average 
down the noise and fit the foregrounds. However, compli- 
cations arise when processing multiple pixels. Different 
pixels come from different lines of sight, so their 21cm sig- 
nals are either slightly different realizations or completely 
independent realizations of equation |(7J, depending on 
how far apart the pixels are from one another. Further- 
more, for pixels that are close to one another, they have 
slightly different signals in general, but on large scales, 
the signals within these pixels are more or less identical, 
while on small scales, the signals are almost independent 
of signals in other pixels. The detailes of this complica- 
tion would probably best be treated via a detailed 3D 
numerical simulation, where thousands of pixels can be 
simulated and the signals from them tested. Since this 
is beyond the scope of this paper, we will simply illus- 
trate the basic effects by two extreme situations, which 
can also be applied to numerically generated signals. We 
will see from the following examples that our method for 
foreground cleaning still works reasonably well. 

4.2.1. Coherent signal approximation 

For closely separated pixels, wc make the crude as- 
sumption that the line of sight 21cm signal in these pix- 
els are identical (same phase and amplitude as in equa- 
tion (0). This approximation will simplify our calcula- 
tion, yet the procedure of doing foreground removal is 
similar to generally incoherent signals as we discuss in 
the next subsection. 

Since the signals are coherent, the total signal for dif- 
ferent pixels is the summation of the same signal and dif- 
ferent foregrounds and noise. We use the same method 
as described in Subsection 14. II to remove the foreground 
from total signal along line of sight for each individual 
pixel. We then average all of them in real space and 
obtain the averaged cleaned signal. Figure El shows the 
results before and after cleaning. The top and middle 
panels are plotted for a single pixel with 200mK noise. 
The noise wiggles fast on top of the foreground and dom- 
inates the signal. It is impossible to tell the difference be- 
tween situations with and without 21cm signal (top and 
middle panels). The bottom panel shows the cleaned sig- 
nal and residual by applying our method and combining 
4000 such pixels. Although both foregrounds and noise 
are at a level orders of magnitude higher than the sig- 
nal, the resulting cleaned signal still captures the main 
features of the "true" signal and the residual is well con- 
trolled. 

This confirms that the foregrounds can still be removed 
effectively when the noise is orders of magnitude higher 
than the signal. Our fitting method does not introduce 
additional contamination to the signal, even when fore- 
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Fig. 5. — Same as Figure ISl but with 200 times larger detector 
noise (tr = 200mK). In the bottom panel, the recovered 21cm signal 
(blue curve) and the residual (red curve) is the result of averaging 
4000 pixels. All other curves in the plot are from a single pixel. 




k [h/Mpc] 



Fig. 6. — Same as Figure El but with 200 times larger detec- 
tor noise (cr = 200mK). The coherent case residual contaminants 
(dot-dashed red curve) are FFT from coherently averaged residual 
in Figure|S] (red curve). The incoherent case residual contaminant 
(red solid, dashed, and dotted curves) are computed from incoher- 
ently averaging 40000 pixels. The red solid, dashed, and dotted 
curves are assuming 1.7MHz, 8.6MHz, and 17.2MHz bandwidth, 
respectively. 



grounds with many different spectral shapes are averaged 
together. 

Figure shows fc-space signal power spectrum com- 



pared with foregrounds and noise power spectra for single 
pixel and residual power spectrum from averaging 4000 
pixels. Before cleaning, for each individual pixel, the 
signal is completely buried under huge noise and fore- 
grounds. After cleaning, the residual (red dot-dashed 
curve) is orders of magnitude below both signal and orig- 
inal contaminants except on very large scales. Similar to 
Figure^! the plot suggests that we need frequency resolu- 
tion substantially smaller than O.lMHz to take advantage 
of foreground removal sensitivity. 

The number of pixels needed to average down the noise 
varies sharply with the actual noise level. To achieve a 
similar level of accuracy shown in Figure for noise 
~ 500mK, we need approximately 20000 to 30000 pixels. 
For noise ^ lOOmK, we need around 1000 pixels. If the 
noise is about the same level as the signal a — 20mK, we 
only need 50 to 100 pixels to adequately lower the noise. 

When larger numbers of pixels getting combined, sev- 
eral side effects appear such as the variation of spec- 
tral indices among different pixels, signal and foreground 
angular correlation, overestimation of the sensitivity at 
small scales, etc. We will discuss some of these issues a 
bit more in the next subsection. The best approach, how- 
ever, is to combine this method with a ngular approach 
and remove the contaminants in 3D fMcQ uinn et al\ 
2005; Morales. Bowman. & Hewitt 2005). Although this 
is beyond the scope of the present paper, we hope to 
address this further in future work. 

4.2.2. Incoherent signals 

For pixels that are far apart, the line of sight 21cm 
signal in these pixels are no longer the same. They have 
different phases and amplitudes and therefore are more 
or less independent to one another. Here we study the 
case when all signals are completely independent. Com- 
pared to the case studied in the previous subsection, this 
case is the opposite extreme. Signals obtained from real 
observations will probably have a behavior intermediate 
between these tow extremes. Simulating such signals nu- 
merically, our cleaning technique is likely to give residual 
contaminant levels that are intermediate between those 
we find here for these two extreme cases. 

When the signals are incoherent, we still use the same 
method to remove foregrounds from the total signal for 
each individual pixel. However, instead of averaging 
them in real space, we FFT the signal in each pixel to 
Fourier space to obtained the cleaned power spectrum for 
each pixel. Then we average all individual power from 
different pixels and get the final average cleaned power 
spectrum. 

Comparing with averaging coherent signals in real 
space, averaging incoherent signals in Fourier space re- 
quires larger numbers of pixels to remove the foregrounds 
effectively. Figure shows true 21cm power spectrum 
compared with residual foregrounds and noise power 
spectra (red solid, dashed, and dotted curves), defined as 
the difference between average cleaned power spectrum 
and true 21cm power spectrum, from incoherently aver- 
aging 40000 pixels. The noise and foregrounds levels are 
kept the same as in previous coherent example. Although 
we average 10 times more pixels here than for coherent 
case, the residual contaminant is at a level higher than 
the previous case. However, the residual is still reason- 
able. On most of the scales, the residual contaminant is 
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less than 10% of the signal. Also notice in this case the 
21cni power spectrum is best measured for scales around 
k = 0.1/i/Mpc, another consequence from incoherence 
averaging. In previous two examples, namely low noise 
and high noise coherence signal examples, we recovered 
signal instead of power spectrum. 

The three different residual curves in the plot are com- 
puted assuming 1.7MHz (red solid curve), 8.6MHz (red 
dashed curve), and 17.2MHz (red dotted curve) band- 
width, respectively. (We used 1.7MHz bandwidth for 
all previous calculations and figures.) As the band- 
width increases, the residual power decreased especially 
on smaller k. That is to say a larger bandwidth will help 
foreground removal at large scales. 

So the bottom line is, for incoherent signals, our 
method for foreground cleaning still works, yet its effi- 
ciency is reduced due to the fact that signals are indepen- 
dent. In this case, we could consider increase bandwidth, 
combine with angular direction measurements, etc., to 
improve efficiency and remove foregrounds effectively. 

5. DISCUSSION 

We have explored how well foreground contamination 
can be removed from a 21 cm tomography data cube by 
using frequency dependence alone. We found that with 
realistic experimental sensitivities, 21 cm tomography is 
limited mainly by foregrounds for scales k <ti 1/i/Mpc 
and limited mainly by noise for k 3> 1/i/Mpc, a result 
which is rather robust to changing the foreground as- 
sumptions. In optimizing the design of upcoming exper- 
iments, a useful rule of thumb is therefore to make the 
channel width substantially smaller than 0.1 MHz, allow- 
ing one to take full advantage of the detector sensitivity 
by pushing residual foregrounds down to the noise level. 
Fortunately, attaining such narrow channel width is re- 
alistic for upcoming experiments, where the analysis is 
all done in dedicated high speed electronics. 

We used a simple "blind" removal technique using no 
prior information about the nature of the foregrounds, 
merely fitting out a quadratic polynomial in log-log for 
the frequency dependence separately for each pixel in 
the sky. The basic reason that this works so well is that 
the foregrounds have a much smoother frequency spectra 
than the 21 signal. 

Although highly effective, this frequency-based clean- 
ing should be viewed as merely one of three comple- 
mentary foreground countermeasures. First, bright point 
sources can be identified as strong positive outliers, and 
the corresponding sky pixels can be discarded since they 
constitute only a small fraction of the total survey area. 
Second, after the frequency-based cleaning step, noise 
and signal can be further distinguished by their differ- 
ent angular correlations, as d escribed in Santos et al. 
ijSantos. Coorav. fc Knoxll2004j) . This angular approach 
will be particularly helpful for early 21 cm experiments 
where the signal-to-noise ratio is limited. The angular 
and frequency-based approaches are therefore comple- 
mentary, and the combination of the two will give the 
best cleaned 21cm signal with which to study the "dark" 
epoch of reionization. 

Although our results are quite encouraging for the 
prospects of doing cosmology with 21 cm tomography, 
much work remains to be done on the foreground prob- 
lem, and we conclude by mentioning a few examples. 



A key assumption is this paper is that the foregrounds 
are dominated by emission mechanisms producing fairly 
smooth spectra. The basis for this assumption is that 
typical atomic and molecular transitions that can pro- 
duce spectral lines correspond to much higher frequencies 
than those relevant to 21 cm tomography. One loophole 
that needs to checked quantitatively is the possible con- 
tribution of recombination lines from hydrogen cascad- 
ing down though very large energy quantum numbers 
n ^ 10^, although early estimates suggest that this is 
not a significant contaminant ('e.o.. .Oh fc Mack .2003}) . 

We performed our analysis on simulated data over a 
small redshift range, limited by our linearization approx- 
imation. It is clearly worthwhile repeating our analysis 
with a proper hydrodynamical simulation of the 21 cm 
signal over the full relevant redshift range. In this case, 
our 3-parameter foreground fit should be generalized to 
one that assumes that the foregrounds are simple only 
locally in log-frequency space. An obvious generaliza- 
tion of our method would be to increase the order of the 
log-log polynomial beyond two. However, we effectively 
want to high-pass filter the observed frequency spectrum 
to clean out foregrounds, and high-order polynominals 
can in principle spoil this by having sharp localized fea- 
tures. A better generalization of our method to long fre- 
quency baselines may therefore be either a cubic spline 
in log-log or a Fourier series expansion. Such end-to- 
end simulations will also be a valuable tool for quanti- 
fying how redshift space distortions (whereby the pecu- 
liar velocity of the gas breaks the one-to-one correspon- 
dence between redshift and frequency) can be exploited 
to separate the effect s of the matter power spe ctrum from 
the "gastrophysics" ffiarka na & Loehll2n04bl) . This be- 
comes important esp ecially for channel width < O.lMHz 
ijDesiacaues fc Nusseili2004t llliev et a/.l |2002b') . 

In summary, the potential science return from 21 cm 
tomography is enormous, both for understanding the 
reionization epoch and for probing inflation and dark 
matter with precision measurements of the small-scale 
power spectrum. Our calculations strengthen the conclu- 
sion that foreground contamination will not be a show- 
stopper. The current situation is similar to the quest for 
the cosmic microwave background in the 1980's in that 
the cosmological signal has not yet been detected, but 
better in the sense that the amplitude of both signals 
and foregrounds are approximately known, guaranteeing 
success if the engineering challenges can be overcome. 
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